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Abstract 

In a frequency selective slow-fading channel in a MIMO system, 
the channel matrix is of the form of a block matrix. This paper pro- 
poses a method to calculate the limit of the eigenvalue distribution 
of block matrices if the size of the blocks tends to infinity. While it 
considers random matrices, it takes an operator-valued free probabil- 
ity approach to achieve this goal. Using this method, one derives a 
system of equations, which can be solved numerically to compute the 
desired eigenvalue distribution. 

The paper initially tackles the problem for square block matrices, then 
extends the solution to rectangular block matrices. Finally, it deals 
with Wishart type block matrices. For two special cases, the results 
of our approach are compared with results from simulations. The first 
scenario investigates the limit eigenvalue distribution of block Toeplitz 
matrices. The second scenario deals with the distribution of Wishart 
type block matrices for a frequency selective slow-fading channel in 
a MIMO system for two different cases of ur = tit and = 2tit- 
Using this method, one may calculate the capacity and the Signal-to- 
Interference-and-Noise Ratio in large MIMO systems. 
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1 Introduction 



With the introduction of some sophisticated communication techniques such 
as CDMA (Code-Division Multiple-Access) and MIMO (Multiple-Input Mul- 
tiple-Output), the communications community has been looking into analyz- 
ing different aspects of these systems, ranging from the channel capacity to 
the structure of the receiver. It has been shown that the channel matrix plays 
a key role in the capacity of the channel PjJ |2| as well as in the structure of 
the optimum receiver [3J Ej. More precisely, the eigenvalue distribution of 
the channel matrix is the factor of interest in different applications. 

For a MIMO wireless system with ny transmitter antenna and ur receiver 



antenna, the received signal at time index n, Y n = [yi <n , 
as follows: 



j ynn,n] 



will be 



Y n = HX n + N n , (1) 

where H is the channel matrix, X n = [xi )Tl , • ■ • , x nTi „] T is the transmitted sig- 
nal at time n and N n is the noise signal. The channel matrix entries hijs re- 
flect the channel effect on the signal transmitted from antenna j in the trans- 
mitter and received at antenna % in the receiver. In a more realistic channel 
modeling, one may consider the Intersymbol-Interference (ISI) Chapter 2]. 
In this case, the channel impulse response between the transmitter antenna j 

h { : j) ■ ■ ■ h ( : j \ h { : j) 1 ' 



and the receiver antenna i is a vector hij 



< l l ""2 '"L-l '"L 

where L is the length of the impulse response of the channel (number of the 
taps). Consequently, the channel matrix for a signal frame of K will be as 
follows: 





' A, 


A 2 




A L 










" 







A 1 


A 2 




A L 










H = 








A, 


A 2 




A L 


■ 






































At 


A 2 ■ 


■ A L _ 



(2) 



sec 



where there are K — 1 zero-matrices in each row and Ai = (hf ) i=\ t ... , nfl 

i=i,- ,n T 

Fig. ^ for the block diagram). To calculate the capacity of such a channel, 
one needs to know the eigenvalue distribution of the H*H 6J. 

Free probability and random matrix theory have proven to be reliable 
tools in tackling similar problems [U]. Tse and Zeitouni [7j applied random 
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matrix theory to study linear multiuser receivers, Moustakas et. al. jH] 
applied it to calculate the capacity of a MIMO channel. Miiller jU] employed 
it in calculating the eigenvalue distribution of a particular fading channel 
and later Debbah and Miiller JU] applied it in MIMO channel modeling. 

While the most basic random matrices are Gaussian and Wishart random 
matrices, and it is now quite widely known (see, for example, jB]) how to use 
tools from free probability (like the R- or S'-transform ^T]) to deal with such 
random matrices, there are many interesting situations, as the matrix H 
above, which cannot be modeled with these simple kind of random matrices. 
Very often these can, at least in first approximation, be assumed to be block 
random matrices. These are matrices built out of blocks, where each block 
is a Gaussian random matrix; however, the variance of the entries might be 
different from block to block and, in particular, some of the blocks might 
repeat at different places in the matrix. Our main example for such block 
matrices is the channel matrix H as in (J2J), but there are also other interesting 
matrices of this kind. In particular, we will also consider block Toeplitz 
matrices, like the 3x3 block matrix 



where A,B,C are independent selfadjoint (which we use as a synonym of 
Hermitian) NxN Gaussian random matrices. Such matrices were considered, 
e.g., in |12j . but their limiting eigenvalue distribution for N — > oo remained 
open. 

While Miiller ^3] has addressed a similar problem using free probability, 
in this paper we will show how to use a more general version of free probabil- 
ity theory, so-called operator-valued free probability theory, to calculate the 
asymptotic eigenvalue distribution of such block matrices. It was the funda- 
mental insight of Shlyakhtenko that operator- valued free probability 
can be used for dealing with special block matrices. Our main theorems 
could, by building on Shlyakhtenko's work, be derived quite directly by com- 
bining various results from the literature on operator-valued free probability 
theory. However, since most readers will not be familiar with operator-valued 
free probability theory, we prefer to give also a more direct proof which does 
not assume any prior knowledge on operator-valued free probability. Actu- 
ally, the "operator- valued" structure of the result will arise quite canonically 
in our proof, and thus our presentation will hopefully convince the reader not 
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only of the power but also of the conceptual beauty of operator-valued free 
probability theory and serve as an appetizer to read more about this field. 

In Section El some basic notations and definitions that are used through 
the paper are introduced JE]- The theorems are stated in Section|3]while they 
are proved in Section 0] In Section |3] it is shown as a warm up exercise how 
the Marchenko-Pastur law can be derived using the proposed method. Then, 
two general scenarios are studied in this section. In the first scenario the 
eigenvalue distribution of block Toeplitz matrices is computed and numerical 
results are shown for 3 x 3, 4 x 4 and 5x5 block Toeplitz matrices. The 
second example computes the eigenvalue distribution for a MIMO system 
with Intersymbol- Interference (ISI) in two different cases. Section|S]concludes 
with a discussion on the convergence and numerical behaviour of the proposed 
method. Conclusion remarks are drawn in Section [HJ 

2 Preliminaries 

2.1 Notations 



The following notations are adopted in the paper: 





Tensor (Kronecker) product of matrices 




direct sum 


3(1) 


Imaginary part of X 


h 


d x d Identity matrix 


X 


complex conjugate of X 


5ij 


Dirac delta function 


X* 


Hermitian conjugate of matrix X 



2.2 Gaussian family and Wick formula 

The entries of our random matrices will consist of Gaussian random variables; 
since entries might repeat at various places in the matrix, not all entries 
are independent and we use the notion of a Gaussian family to describe 
this general situation. A Gaussian family is a family of random variables 
X\,...,x n which are linear combinations of independent Gaussian random 
variables. Clearly, the whole information about the joint distribution of such 
a family is contained in the covariance matrix C = (%)^ J=1 , E [xiXj] = Cij. 
There is a very precise combinatorial formula, usually called Wick formula, 
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which allows to express the higher joint moments in a Gaussian family in 
terms of second moments (i.e., in terms of the covariance matrix). Namely, 
for all k G N and 1 < i (l) , • • • , i (k) < n, we have: 

E [xni) • • ■ Zi(fc)] = ^2 II E [x i(r )X i(s) ] , (3) 

where Vz(k) is the set of all pairings of the set {1, • • • ,k}. Since there 
is no pairing of an odd number of elements, this formula also contains the 
statement that all odd moments of a Gaussian family vanish. 



2.3 Cauchy transform and Stieltjes inversion formula 

Our main analytical object for determining the eigenvalue distribution of a 
random matrix is the Cauchy transform of a probability measure fi on R. 
It is defined as follows: 

GM= l —Jn{z) t zeC + , (4) 

JR z ~ 1 

where C + := {s + it\s,t G R, t > 0}. is analytic and takes values in 
C~ := {s + it\s,t G R, t < 0}. For a compactly supported /i, G M can be 
expanded as a power series: 

oo 
n=0 

where r := sup {|i| | t G supp (/i)} and a n := f R t n d[i (t) is the nth moment 
of \x for n > 0. Using (jij), it is clear that the Cauchy transform has the 
following property: 

lim zG^ (z) = 1. (6) 

2;GC+,|2;|^oo 

One can recover the probability measure \x from its Cauchy transform by the 
following formula, known as the Stieltjes inversion formula: 



where 



da (t) = lim h e (t) dt 



h e (t) := --SGJt + ie) , VtGR. (7) 

71 
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3 Statement of the Theorems 



3.1 Selfadjoint block matrix with square blocks 

We will consider d x rf-matrices with a given block structure; the blocks 
themselves will be large Gaussian N x iV-random matrices; our aim is to 
calculate the limit eigenvalue distribution for such block matrices if the size 
iV of the Gaussian matrices tends to infinity. The main point is that the 
blocks are not all different, but they might repeat (as themselves or their 
adjoint) at different places. Otherwise from that they should be independent. 
We write our matrices as 

d 

x N = Y. E n® A{i,j) i ( 8 ) 

where Eij are the elementary d x <i-matrices having a f at the entry (i, j) and 
otherwise. For each i,j — l,...,dwe have a Gaussian N x N random matrix 
A^' j \ which we call a block. In each block all the entries are Gaussian 
and independent. Furthermore, we want to be selfadjoint, which means 
that A( l 'fi = A^'^* for all i,j = l,d. In particular, the blocks A^ on the 
diagonal must be selfadjoint. We would like to leave it open whether the 
blocks A^ with % ^ j are selfadjoint or not - this depends on the concrete 
situation. 

A convenient way to encode the information about which blocks agree, is 
by specifying the covariance cr(i,j; k, I) between an entry in the matrix A^'^ 
and an entry in the matrix A^ k ' l \ In the case that the two blocks are either 
independent or the same (where we also include the possibility that one is 
the adjoint of the other), the covariance a takes just on the values or 1. 
However, we can also be more general and allow arbitrary cr's. The only 
constraint comes from the requirement that Xn is selfadjoint which means 
that we must have 

v(i,j;k,l) = a(k,l;i,j) (9) 

for all i,j,k,l — l,..., d. 

If we denote the entries of the matrix A^^ by with r,p — 1, . . . , N, 
then we can summarize all the above by requiring that the collection of all 
entries {aip^ \ i,j = 1, . . . , d, r,p = 1, . . . , N} of the matrix X N forms a 
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Gaussian family which is determined by 

= a p j r^ for alH, j = l,...,d, r,p = 1,...,N 
and the prescription of mean zero and covariance 

E[af^a^] = -8 rs 5„-a(i,j;k,l), (10) 

where 

n = dN. 

Note that the above gives also the covariance between a's and a's as 

and thus the joint distribution of the matrix entries of A is uniquely deter- 
mined. 

Furthermore, with this description we have the possibility of including 
both cases where non diagonal blocks are selfadjoint or not (or mixtures of 
this). Consider a situation where a takes on only the values or 1. Then a 
non-diagonal A^ is selfadjoint if and only if cr(i,j; = 1. As an example, 
contrast the situation 

"0 A x " 

A t Ai 
A x 

where A 1 is a selfadjoint Gaussian iV x iV-matrix with the situation 

"0 B x " 

B\ B 1 , 
B{ 

where i?i is a non-self adjoint Gaussian iV x iV-matrix. 

In the first case a(i,j;k,l) = 1 for all 16 possible combinations of 
and (k, I) from the set {(1, 2), (2, 1), (2, 3), (3, 1)}, whereas in the second case 
<r(i,j] k, I) is only non-vanishing for 

<r(l, 2; 2, 1), <r(l, 2; 3, 2), a(2, 3; 2, 1)^(2, 3; 3, 2) 

and the "adjoint values" 

a(2, 1; 1, 2), a(3, 2; 1, 2), <r(2, 1; 2, 3), <r(3, 2; 2, 3). 
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Notation 1. Fix a natural number d and a covariance function 

such that © holds. 

1) Md(C) are just the d x <i-matrices with complex entries, 

M d (C) : {(d,^ , /../ 1 d\. 

2) We define a covariance mapping 

V : M„(C) -> M d (C) 

as follows: 

For D = (dy)&=i G M d (C) we have Typ) = (7 7 (Z?) ii )^ =1 with 



:= d E cr(i,k;l,j)d kl . 

k,l=l 

3) Furthermore, tr^ denotes the normalized trace on M^(C), i.e., 

1 d 

tr <^) : j !>>]«■ 

Theorem 2. Wii/i i/ie above notation, for each JVgN, consider block ma- 
trices (JEJ), where, for each i,j = l,...,d, the blocks = {o>rp^) rp=1 
are Gaussian N x N random matrices such that the collection of all entries 
{aip^ | i, j = 1, . . . , d, r,p — 1, . . . , iV} of the matrix X N forms a Gaussian 
family which is determined by 



afyfi = dpr for alli,j = 1, . . . , d, r,p = 1, . . . , N 

and the prescription of mean zero and covariance (ll(Jj) . 

Then, for N — > oo, the n x n matrix has a limiting eigenvalue distri- 
bution whose Cauchy transform G(z) is determined by 

G(z)=tr d (g{z)), 

where Q{z) is an M^(C) -valued analytic function on the upper complex half 
plane, which is uniquely determined by the facts that 

lim zQ(z) = I d , (11) 

|z|-»oo,3f(«)>0 
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and that it satisfies for all z in the upper complex half plane the matrix 
equation 

zG{z)=I d + r 1 {G{z))-G{z). (12) 

The proof of this theorem is just a few lines if one cites the relevant 
literature from operator-valued free probability theory. We will indicate this 
below. However, since most readers will not be familiar with operator-valued 
free probability theory, we prefer to give in the next section a more direct 
proof from scratch. As one will see, the "operator- valued" structure of the 
result will arise quite canonically and the next section is also intended to 
introduce the reader to some relevant notions and ideas of operator-valued 
free probability theory. 

Let us now give, for the reader familiar with operator-value free proba- 
bility, the condensed version of the proof of our Theorem El First, one has 
to observe that the joint moments of the blocks converges to a semi-circular 
family, thus the wanted limit distribution of is the same as the one of a 
d x d-matrix S, where the entries of S are from a semi-circular family, with 
covariance a. By using the description of operator- valued cumulants of this 
matrix in terms of the cumulants of the entries of the matrix (see [T7| ), it is 
obvious that S is a M^(C)-valued semi-circular element, with covariance 77. 
The equation for G{z) follows then from the basic /^-transform or cumulant 
theory of operator- valued free probability theory, see (HI EH] • 

In the special case that blocks of our block matrix Xn do not repeat (apart 
from the symmetry condition Xn = X^) one might call the block matrix a 
band matrix (in such a situation one often lets d also go to infinity). The 
fundamental observation that limits of Gaussian band matrices are operator- 
valued semicircular elements - and thus that operator-valued free probability 
can be used for determining their eigenvalue distribution - was made by 
Shlyakhtenko in ^3] . There he proved actually the special case of the above 
Theorem |21 for band matrices. In this case, the covariance function 77 maps 
diagonal matrices to diagonal ones and one can consider the block matrix as 
a semi-circular element over the diagonal matrices. 

3.2 Selfadjoint block matrix with rectangular blocks 

In many applications one is also interested in situations where the blocks 
themselves might not be square matrices, but more general rectangular ma- 
trices. Of course, the sizes of the blocks must fit together to make up a big 
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square matrix. This means that we replace n = dN by a decomposition 

n = Nx + 1- N d , and the block will then be a iVj X A^-matrix. We 

are interested in the limit that Ni/n converges to some number ctj. It turns 
out that one can modify our Theorem |21 (and also its proof) quite easily from 
the square to the general rectangular situation. Let us first introduce the 
generalizations of our relevant notations from the square case. Note that 
dependent rectangular blocks can be re-cut into different nonequivalent con- 
figurations of dependent blocks. We will assume that such repartitioning has 
already been done and resulted in the covariance function a(i,j;k,l) that 
can only be different from zero if the size of the block A^'^ fits (at least in 
the limit n — > oo) with the size of the block A^ k ' l \ 

Notation 3. Fix a natural number d and a <i-tuple a = (a\, . . . , a d ) with 
< aij < 1 for all * = 1, . . . , d and a\ + • • • + ad = 1. Furthermore, let a 
covariance function o = (cr(i, j; k, I)). k z _ 1 be given with the property that 
(0) holds and in addition k, I) = unless ai = ai and <x,- = a^- Then 

we use the following notations. 

1) M a (C) are those matrices from M^(C) which correspond to square 



blocks, 



M a (C) := {(dij) G M d (C) | d i:j = unless = aj}. 
2) We define the weighted covariance mapping 



Va : M a (C) -> M a (C) 



as follows: 



For D = (dtj) 




G M a (C) we have T] a (D) = (r)a(D)ij) 




with 



d 




k,l=l 



3) Furthermore, the weighted trace 



tr a : M a (C) -> M a (C) 



is given by 



d 




i=i 
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Theorem 4. With the above notation, for {Ni, . . . , Nj} C N consider block 
matrices 

d 

X Nl _ Nd =Y,E^®A^\ 
*,i=i 

Here Eij are the elementary dxd-matrices having a 1 at the entry (i,j) and 
otherwise and, for each i,j = 1, . . . , d, the A^'i> are Gaussian NiXNj random 

matrices, (o!rp^ r=i,...,H i ■ The latter are such that the collection of 

p =i,..., Nj 

all entries {arp \ i, j = 1, . . . , d, r = 1, . . . , N,p = 1, . . . , Nj} of the matrix 
XN u ...,N d forms a Gaussian family which is determined by 



= afr for alii, j = 1, . . . , d, r = 1, . . . , Ni, p = l,...,Nj 
and the prescription of mean zero and covariance 

E[a^ j) a^ l) ] = ^5 rs 5 pq ■ a{i,j; k, I), 

where we put 

n:=N l + ---N d . 

Then, for n — > oo such that 

lim — - = cti for all i = 1, . . . , d, 



the matrix X Nl ^^ jNd has a limiting eigenvalue distribution whose Cauchy 
transform G(z) is determined by 

G(z)=tv a (g(z)), 

where Q(z) is an M a (C) -valued analytic function on the upper complex half 
plane, which is uniquely determined by the facts that holds and that it 
satisfies for all z in the upper complex half plane the matrix equation 

zG{z)=I d + r la {G{z))-G{z). (13) 

In principle, the statement for the rectangular situation can be reduced 
to the case of square blocks by cutting the rectangular blocks into smaller 
square blocks (at least asymptotically); thus Theorem HI can be deduced 
from Theorem |21 However, we find it more instructive to adapt our proof 
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of Theorem |2] to the rectangular situation - this makes the structure of the 
weighted covariance mapping and the appearance of the weighted trace more 
transparent. We will present this proof of Theorem 0] also in the next section. 
Let us also mention that the considerations of Shlyakhtenko for Gaussian 
band matrices were extended to rectangular matrices by Benaych-Georges 
[2Tij : thus some of his results are closely related to the special case of band 
matrices in our Theorem 0] 



3.3 Wishart type block matrices 



Another type of block matrices appearing in applications is of a Wishart type. 
Here the block matrix H itself is rectangular and one is interested in the 
eigenvalue distribution of HH* . Let us write H = (A^)i=i,...,r, where each 

j = l,...,s 

block A^ ,: >> is a Gaussian Mi x Nj random matrix. Put M := M 1 + ■ — h M r 
and N : = Ni + ■ ■ ■ N s , so H is an M x iV-matrix. The calculation of the 
eigenvalue distribution of HH* can be reduced to the situation treated in 
the previous section by the following trick. Consider 



X 



H 

H* 



(14) 



With n = M + N, this is a selfadjoint n x n-matrix and can be viewed as 
a d x d-block matrix of the form considered in the previous section, where 
n = M + N is decomposed into d := r + s parts according to 

M + N = Mi + • • • + M r + JVi + • • • N s . 

Thus we can use Theorem 0]to get the eigenvalue distribution of X if the Mi 
and Nj go to infinity with a limit for their ratio. 

The only remaining question is how to relate the eigenvalues of X with 
those of HH*. This is actually quite simple, we only have to note that all 
the odd moments of X are zero and 



X 2 



HH* 






H*H 



Thus the eigenvalues of X 2 are the eigenvalues of HH* together with the 
eigenvalues of H*H. Furthermore, note that HH* is an M x M and H*H 
is an N x N matrix. Assuming that M < N (otherwise exchange the role of 
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H and H*) we have then that the eigenvalues of H*H are the eigenvalues of 
HH* plus N — M additional zeros. In terms of the unnormalized trace Tr 
this just means 

Tr(X 2fe ) = Ti((HH*) k ) + Ti((H*H) k ) = 2Tr((M*) fc ) 

Rewritten in normalized traces this gives 

tr n (X 2k ) = — tr M ((HH*) k ) ■ 

n 

Putting all these together, it results in the following relation between the 
Cauchy transform Gx? of X 2 and the Cauchy transform Ghh* of HH*: 

r ( \ M + N n f \ N-M l 

Ghh*\z) = — — Gx^(z) — . 15 

v ; 2M y ' 2M z K J 

Finally, we should also rewrite our equation for the Cauchy transform Gx of 
X in terms of Gx 2 - Since X is even, both are related by 

z-G x i(z 2 ) = Gx(z). 

By noting that the operator- valued version Q(z) depends only on z 2 , we can 
introduce a quantity 7i by 

z-H(z 2 ) =G(z). 

Then we have 

lim G x *(z) =tx a [H(z)}, (16) 

n— >oo 

and the equation (JT2J) for Q becomes 

zH(z)=I d + zr ] (H(z))-H(z). (17) 

The equations (fTH|) . (fTTj) . together with (fTKj) determine the asymptotic eigen- 
value distribution of HH* by 

lim G HH *(z) = 6ti a H(z) - — , (18) 

M+N->oc Z 

where 

ft - ]irn M + N - 1 ft - ,;„, N Z M - * ~ <** 
u = lim = — =t , on = lim = =7: . 

n^oo 2M 2jJ j=1 aj tv^oo 2M ^T>j=i a i 
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Note that the subtraction of a pole at in (|15|) reflects just the occurrence 
of spurious zeros in X 2 which are coming from H*H. 

It would seem more adequate to have access to the information about the 
eigenvalues of HH* without having to invoke the eigenvalue distribution of 
H*H. To simplify notations we will restrict in the following to the situation 
where all blocks A^' are square matrices i.e., Ni — • ■ ■ — N r — M 1 = ■ ■ ■ = 
M s , and thus N = rN± and M — sN^ (Note that H itself might still be 
rectangular, i.e., we are allowing that r ^ s.) One can of course write our 
(r + s) x (r + s) matrix 7i 2-block matrix 



H(z) 



Qi(z) g 3 (z) 



where Qi(z), Q 2 (z), and G^z) are rxr-, sxs-, and r x s-matrices, respectively. 
From our considerations in the next section is fairly easy to see that 



lim G H h*{z) = tr r (£i(»). 



(19) 



One should note that the quadratic matrix equation (|T7|) for 7i will in general 
couple all the various Qi, so that even in this formulation the calculation of 
the distribution of HH* will still involve H*H. However, there are special 
situations where one can eliminate H*H quite easily from our equations. 
We will examine one such situation, which appears in most applications to 
wireless communications, in more detail in the rest of this section. 

Let us consider the special version of our Wishart type block matrices 
where the blocks of H themselves are non-selfadjoint Gaussian random ma- 
trices (which are still, as before, square matrices). In this case the covariance 
a for the blocks of ()14j) couples effectively only blocks from H with blocks 
from H* and thus r] : M r+s (C) -> M r+d (C) maps M r (C) ©M S (C) to itself by 



77 : 



~D l " 




m(D 2 ) 





D 2 










where 



77! : M r (C) -> M,(C) and r] 2 : M S (C) M r 
One sees now easily that H must take on values in M r (C) © M S (C), thus 



H(z) 









G 2 {z) 
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where Qi and Q 2 are M r (C) -valued and M s (C)-valued, respectively, analytic 
functions in the upper complex half plane. The equation (jT7j) for Tt splits 
now into the two equations 

zQi(z) = I r + zt) 2 (G 2 (z)) ■ Qx(z) 

and 

zQ 2 {z) = I s +ZT)i(Gi(z)) • Qliz). 

One can eliminate Q 2 from those equations by solving the second equation 
for Q 2 and inserting this into the first equation, yielding 

Z g 1 (z) = i r + m({is - vi(Gi{z))Y l ) ■ Gib) 

Together with f)19|) this gives directly the eigenvalue distribution of HH* . 



4 Proof of the Main Theorems 
4.1 Proof of Theorem [2] 

We will prove Theorem 2 by calculating the moments of the averaged eigen- 
value distribution of the n x n random matrices in the limit N — > oo. 
Note that the m-th moment of the averaged eigenvalue distribution is just 
given by the expectation of the normalized trace tr n of the m-th power of 
X]\f. For this we can calculate 

1 d 

E[ti n (X%)] = - ^[trjv(A (i(1) ' i(2)) • • ■ At'W-'t 1 )))] 

,...,i(m)=l 

1 d N 

_ 1 P r (i(l),i(2)) (i(m),i(l))l 

~~ dN ^ L a r(l)r(2) • ' ' a r(m)r(l) J 
i(l),...,i(m)=l r(l),...,r(rri)=l 

Now we need to invoke the Wick formula for calculating the expectation of a 
product of entries from our matrix. The Wick formula allows to express all 
higher moments in random variables from a Gaussian family in terms of the 
covariances by a nice combinatorial formula, by summing over all pairings of 
the appearing factors. In our case this gives 

w v (i(l),<(2)) _(t(ro),t(l))-i _ TT Eir_(<(p),t(iH-l))_(<(«),i(«+l))i 

rj [ a r(l)r{2) "' tt r(m)r(l) J - 11 L«r(p)r(p+1) tt r(q)r( g +l) J" 

7rEV2(jn) (p,q)S7r 
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Here V%[ni) denotes the pairings of a set of to elements. A pairing is just a 
decomposition into subsets (called blocks of the pairing) with two elements. 
The product in the above Wick formula runs over all blocks of the pairing 
7i. If we now use the formula (jlUj) for the covariances of the entries of our 
random matrix Xjy, then we can continue with 

p r (i(l),<(2)) (i(m),i(l))l 
^ L U r(l)r(2) ' ' ' U r(m)r(l) J 

= Yl II n ' °"(*(p)>*(P + 1 )'iK<l)> i ( ( l+ !)) • S r(p)r(q+1) " <*r(p+l)r(«)- 

In order to understand better the contribution of a pairing 7r according to 
this formula, it will be useful to interpret this formula in the language of 
permutations. By S m we will denote the group of permutations on to ele- 
ments. Let us denote by 7 G S m the cyclic permutation 7 = (1,2, ... ,m). If 
we identify a partition n G V-iim) with a permutation in - by declaring 
the blocks to cycles, i.e., 7r(p) = g and 7r(g) = p for (p,q) G 7r - then one 
recognizes easily that 

N N m 

II S r(p)r(q+1) ■ 8r(p+l)r(q) = ^ II 6 r(p)r(yn(p)) = N* 1 ^ , 
r(l),...,r(m)=l (p,<j)G7r r(l),...,r(m)=lp=l 

where #7?r denotes the number of cycles of the permutation 771-. 

With this we can continue our above calculation of the expected trace of 
X™ as follows: 

E[\x n (X$)\= n *^~ m/2 ~ 1 ^ E II CT(i(p),i(p+l);i(q),i(q+l 

n£V2(m) i(l),...,i(m)=l (p,<j)e7r 

Note that for d = 1 (and cr(l, 1; 1, 1) = 1) we recover the well-known 
genus expansion 

E[ti n {A m )} = J2 n*^- mt2 - x 
of a selfadjoint Gaussian n x n random matrix. 



It is quite easy to see that 

#77r — m/2 — 1 < for any 7r G Vi(wn) 
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and that we have equality there exactly for the so-called non-crossing pairings 
7i G NC2(m). The latter are those pairings tt of 1, 2, . . . , 2m for which the 
blocks of tt can be drawn in such a way below the numbers that they do not 
cross. Another useful way of describing non-crossing pairings is the following 
recursive procedure: A pairing of 2m numbers is non-crossing if and only if 
one can find a block consisting of two consecutive numbers such that after 
removing this block one gets a non-crossing pairing of the 2m — 2 remaining 
numbers. For illustration, here are the 5 of the 15 pairings of 6 elements 
which are non-crossing: 



1 2 3 4 5 6 

u u u 



1 



u 



2 3 4 5 6 

u 



1 



2 3 4 5 6 

u 



u 



and 



So in the limit, N 
formula and we get 



2 3 4 5 6 

u u 



1 2 3 4 5 6 

u 



oo, only the non-crossing pairings survive in our 



lim E[ti n (X%)] 



7r6WC2(m) 



where 



: d m/2+l 



Yl II < 7 ( i (p)>*(p + 1 ); i (9)> i (9 + 1 )) 

i(l),...,i(m)=l (p,q)e7r 



In the case d — 1, all K,^ would be equal and we would recover the 
moments of a semicircular distribution (in the form of the Catalan num- 
bers counting non-crossing pairings). This would reproduce more or less the 
derivation of Wigner of his famous semicircle law. For d > 1, however, the 
)C n are in general different for different tt, and, even worse, there does not 
exist a straightforward recursive way of expressing JC n in terms of " smaller" 
K, c . Thus we are outside the realm of the usual techniques of free probability 
theory. 

However, one can save most of those techniques by going over to an 
"operator- valued" level. The main point of such an operator- valued approach 
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is to write JC n as the trace of a d x d- matrix K n , and then realize that k w has 
the usual nice recursive structure of semicircular elements. 
Namely, let us define a matrix 

^7T = j))ij = l 

by 

d l 

K n(hj) '■= Yl *»(i)^-i(m+i) IJ 2 •a(i{p),i{p+ l);i(g) ■ 

j(l)...,i(m),i(m+l)=l (p,<?)e7r 

Then clearly we have 

= tr^). 

Furthermore, the value of be determined by an iterated application 

of our covariance mapping rj : M d (C) — > M rf (C). Namely, the value of K n is 
given by an iterated application of this mapping r\ according to the nesting 
of the blocks of ir. If one identifies a non-crossing pairing with a putting of 
brackets, then the way that rj has to be iterated is quite obvious. Let us 
clarify these remarks with an example. 
Consider the non-crossing pairing 



u 



u 



7r = {(l,4),(2,3),(5,6)}e7VC 2 (6) = 
The corresponding 

l d 

M*, j) = # E °^ ? ( 2 ); *( 4 )> *( 5 )) • o-(i(2),i(3); i(3), i(4)) 



(i 3 

i(2),i(3),i(4),i(5),i(6)=l 



.o-(z(5),i(6);z(6),j). 



7r is a non-crossing pairing of the numbers 1, 2, 3, 4, 5, 6. Let us put the 
indices % = i(l), i(2), i(3), i(4), i(5), i(6), i(7) = j between those numbers. 
Each block of it corresponds to one factor a in K n and the arguments in this 
a are the indices to the left and to the right of the two numbers which are 
paired by this block. 

. 1 i(2) 2 i(3) f i(4) 4 i(5) 5 i(6) 6 j 



18 



As mentioned before, for a non-crossing pairing one always has at least one 
block consisting of neighbours; in our case such a block is (2, 3) G n (actually 
(5,6) is another one). This means that we can sum over the corresponding 
index i(3) without interfering with other blocks, giving 



K 7r(2) J J 



1 d 



a(i,i(2);i(4),i(5)) ■ a(i(5),i(6);i(6),j) 



i(2),i(4),i(5),i(6)=l 



- d a(<(2), i(3);<(3),<(4)) 



i(3)=l 



1 ^ 

= j 2 a ^ l{2 ^ ^( 4 )^( 5 )) • cr(»(5),i(6); i(6),j) • 77(l)i(2)<(4) 



i(2),i(4),i(5),i(6)=l 



Effectively we have removed the block (2, 3) and replaced it by the matrix 
77 (id); we are left with 



i(2) i(4) 

Wd)} i(2)i{4) 



i(5) 



i(6) 



Now the block (1,4) of 7r has become a block consisting of neighbours 
and we can do the summation over i{2) and i(4) without interfering with the 
other blocks, thus yielding 



1 d 

K n (i,j) = - <7(*(5),i(6);i(6),j) 



i(5),i(6)=l 



(2)i(4) 



i(2),i(4)=l 



i £ ^(5),^(6)^(6),j)-^(^(/ d ))] 



ii(5) 



i(5),i(6)=l 



We have now removed the block (1,4) and the effect of this was that we had 
to apply 7] to whatever was embraced by this block (in our case, i](I d )). We 
remain with 



[v(v(i d ))] 



i(5) 



m(5) 



i(6) 



6 



J 
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Now we can do the summation over i(5) and i(6) corresponding to the 
last block (5, 6) of n, which results in 

**(ij)= E [v(v(h))] u{5) (l E a(i(5),i(6);i(6),j)) 

i(5)=l i(6)=l 
cZ 

= E [ ? ?W / «'))]ii(5)^( 1 k5W 
i(5)=l 

= [77(77(7,)) .77(7,)].. 

Thus we finally have 

which corresponds to the bracket expression 

(X(XX)X)(XX). 

In the same way every non-crossing pairing results in an iterated application 
of the mapping 77. For the five non-crossing pairings of six elements one gets 
the following results: 



1 2 3 4 5 6 



U U U |j 



1 2 3 4 5 6 

u 



77(7,) • 77(7,) • 77(7,) 77(7,) -77(77(7,)) 



1 2 3 4 5 6 



u 



u 



and 



1 2 3 4 5 6 

u u 



20 



1 2 3 4 5 6 

u 



Thus for m = 6 we get 
lim E[tr(A 6 )] 

n— >oo 

= tr d |77(J d ) -^(/d) -77( J rf )+77(J d ) -77 (r7(7 d )) +77 (77^)) -?7(^)+?7(^(^) -^(^)) (?? (?7(^))) } - 
Let us summarize our calculations for general moments: We have 
lim E[ti n (X%)] = tr d { V 

7r£NC2(m) 

where d x d matrices, determined in a recursive way as above, given 

by an iterated application of the mapping rj. 

Let us denote by £ the limiting operator- valued moments of X N , i.e., 

S(X m ) := Yl K ^ 

ir£NC2{m) 

so that we have 

lim E[tr n (X%)}=tr d (£(X m )). 

An element X whose operator- valued moments S(X m ) are calculated in 
such a way is called an operator-valued semicircular element. 

The above description of the moments S(X m ) can, similarly as for an 
ordinary semicircular variable, be reformulated in a recursive way which leads 
directly to our equation (JT5J). To get this recursive description, one has to 
observe that if 7r is a non-crossing pairing of m elements, and (1, r) is the block 
of 7r containing 1, then the remaining blocks of 7r must fall into two classes, 
those making up a non-crossing pairing of the numbers 2, 3, . . . , r — 1 and 
those making up a non-crossing pairing of the numbers r + 1, r + 2, . . . , m. 
Let us call the former pairing 7Ti and the latter 7T2, so that we can write 
7r = (1, r) U 7Ti U 7r 2 . Then the above description of k n shows that 
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This results for the operator valued moments in the recurrence relation 

m-2 

£[x m ] = v{£\x k }) ■ £\x m - k - 2 }. 

If we go over to the corresponding generating power series, 

oo 

M{z) = ^2S[X m ]z m , 

m=0 

then this yields the relation 

M(z) = I d + z 2 r](M(z)) ■ M(z). 

Note that M(z) := tid{Ai (z)) is the generating power series of the moments 
E\tT]y(X™)] for N — > oo, thus it is preferable to go over from A4(z) to the 
corresponding operator-valued Cauchy transform 

g( z ) ;= -M(l/z). 

For this the above equation takes on the form 

zG(z) = I d + r ] (g(z))-g(z), 

which is exactly our equation ()12j) . Furthermore, we have for the Cauchy 
transform G of the limit eigenvalue distribution of our block matrix that 

„ t , M(l/z) Mil/z). 

G(z) = y ; = tr d ( y ; ) = tv d (g(z)). 

Let us finally make a few remarks about the support of the limiting 
eigenvalue distribution and the validity of equation (JI2J) for all z in the upper 
half plane. 

First note that the number of non-crossing pairings of 2k elements is given 
by the Catalan numbers which behave asymptotically for large k like 4 fc . (Of 
course, for m odd there are no pairings of m elements at all.) This implies 
that we can estimate the (operator) norm of the matrix S(X m ) by 

\\£(X m )\\ < \\rj\\ m ■ #iVC 2 (m) < c- IMP -2 m , 
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where c is a constant independent of m. Applying tr d , this yields that the 
support of the limiting eigenvalue distribution of is contained in the 
interval [— 2||?y||, +2||r/||]. Clearly, since all odd moments are zero, the mea- 
sure is symmetric. Furthermore, the above estimate on the operator-valued 
moments £{X m ) shows that 

k=0 

is a power series expansion in 1/z of Q(z), which converges in a neighborhood 
of oo. Since the mapping 

D » -I d + -rj(D) ■ D 

z z 

is a contraction for \z\ sufficiently large, Q{z) is, for large z, uniquely deter- 
mined as the solution of the equation (|12p. 

In order to realize that Q(z) can actually be extended to an analytic 
matrix-valued function on the whole upper complex half plane one has to 
note that the equation 

lim E[ti n (X%)] 

n— >oo 

1 d 

Tr£NC2(m) i(l)....,i(m)=l (p,g)S7r 

can also be read in the way that the limiting moments of X^ are given by 
the moments of a d x cf-matrix 

where the entries of X form a circular/semicircular family living in some non- 
commutative probability space (A, <f), with covariance (p(cijCki) = cr(i,j; k, I): 

lim E[ti n (X%)] = V ®ti d {X m ). 

n—*oc 

What we showed above was then essentially that X is an operator-valued 
semicircular variable over M^(C), where the expectation 

£ : M d {A) = M d (C) ®A^ M d (C) 



23 



is given by applying ip entrywise, i.e., 8 = 1 <S> tp. In this language, 



which shows that it is actually an analytic function on the whole upper 
complex half plane. The validity of (fT2*j) for all z in the upper half plane 
follows then by analytic continuation. 



4.2 Proof of Theorem |U 

The proof of Theorem 0] is very similar to the one of Theorem El We will 
concentrate here mainly on the modifications compared to the latter one. 

Again, one calculates the averaged traces of powers by invoking the Wick 
formula. This yields 

d N i(l) N i(m) 

w % Nj)] = i £ E-E E II s 

i(l),...,j(m)=l r(l)=l r(m)=l 7r£"P2 (m) 
• ^(*(p),i(P + l);*(?),*(?+l)) • Sr(p)r(q+1) ' <*r(p+l)r(fl) 
d N 

= E ^~ m/2 - 1 E n^r II ^(<(p),<(p+ i);<(<?).<(<?+ 1))- 

Tr&V2{m) i(l),...,i(m)=l cEyn (p,<?)6t 

Here iVj( c ) denotes the value A^( fc ) which belongs to the cycle c of •yn. Note 
that the appearance of the cx-factors ensures that, although i(k) might be 
different for different k in the cycle c, all iV^) for k G c are actually the 
same. 

Thus in the limit n — > oo, again only the non-crossing pairings survive 
and we get as before 

lim E[tr(X% N )] = E ^ 

weNC2(m) 

where now the contribution of a tt G NC2(m) is given by 

^ := e n n ° r ( i (^)' i (p+ !);%)>%+ 1))- 

i(l),...,i(m)=l c677r (p>9)G7r 
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As before we can write this as a trace of a matrix-valued object - but this 
time everything is weighted with the vector a. Namely, let us define the 
matrix 

by 

d 

i(l)...,i(m),i(m+l)=l 

n a i{c ) n o-(«(p),«(p +!);«(?),*(?+ 1)), 

where this time we only take the product over those cycles of 77r which do not 
contain the element 1. The factor a^i) corresponding to the cycle containing 
1 will be used as a weighting factor for the remaining trace: 

JC n = tr a (K 7r ). 

One can now check easily that the value of is calculated in the same 
way as before, by an iterated application of the (now weighted) covariance 
mapping 7] a , nested according to the nesting of the blocks of the non-crossing 
pairing ir. 

Let us summarize our calculations: We have 

Km E[tr n (X% i _ Nd )]=tr a { ^ ^}, 

TreNC2(m) 

where d x d matrices in M^(C), determined in a recursive way as 

above, and given by an iterated application of the mapping rj a . 
The derivation of the equation (fT3|) is then as before. 

5 Results and Discussion 

Our theorems give us the Cauchy transform G of the asymptotic eigenvalue 
distribution of the considered block matrices in the form G(z) = tr a (Q(z)), 
where Q{z) is a solution to the matrix equation (JT2J), (fT3j). or (JT7J). We 
recover the corresponding eigenvalue distribution \x from G in the usual way, 
by invoking Stieltjes inversion formula 

dfi(x) = lim QG(x + ie)dx, (20) 

7T e\0 
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where the limit is weak convergence of measures (as we will remark in 

sec, inn ■ 

Usually, there is no explicit solution for our matrix equation, so that we 
have to rely on numerical methods for solving those. Note that we do not 
get directly an equation for G. We first have to solve the matrix equation, 
then take the (weighted) trace of this solution. Thus, in terms of the entries 
of our matrix Q, we face a system of quadratic equations which we solve 
numerically using Newton's algorithm |21j . 

In many concrete cases the matrix Q exhibits some special structure (i.e., 
some of its entries might be zero or some of them agree); taking this into 
account from the beginning reduces the complexity of the considered system 
of equations considerably. This special pattern of the entries of Q depends of 
course on the considered block structure of our problem, i.e. on the covariance 
function rj. Let us point out how one can determine this pattern of Q, relying 
on the knowledge of rj. Let us concentrate on the equation (fT^j) for square 
matrices. As we remarked in the previous section, for large z, one can produce 
a solution to this equation by iterating the map 

D h-> -Id + -rj(D) ■ D. 

z z 

Since a pattern in the entries of Q which is valid for large z must, by analytic 
continuation, remain valid everywhere, we can use the above map to deter- 
mine this pattern. Since the factor 1/z plays no role for recognizing such a 
pattern, we start, for example, with the identity matrix D = Id and iterate 
the mapping D i— ► Id + rj(D) ■ D a few times until the pattern in the entries 
of D stabilizes. This is then the wanted pattern for the entries of Q. 

In this section we will consider concrete examples of block matrices and 
compare our solution from Theorem |21 or 0] with the result of simulations. 
First, the well-known Marchenko-Pastur law is derived analytically. Then the 
eigenvalue distribution of two different block matrix scenarios is calculated. 
We conclude this section with remarks on convergence and the numerics of 
the proposed method. 

5.1 Marchenko-Pastur law 

This is a warm-up example in which we re-derive the celebrated Marchenko- 
Pastur law from Theorem 0] 
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5.1.1 Proposition: Marchenko-Pastur 



If Hn is and M x N matrix of i.i.d. random variables of unit variance, then 
as N —>■ oo, N/M — > A > 0, the empirical law of the eigenvalues of HH*/M 
converges weakly with probability one to 



I J(x - 1 - A) 2 -4A 
H{dx) = (1 - Xr5 + — 1 ( i-Va)2< x <(i + Va)2^ ( 21 ) 

Proof. Since this is a well known result, we just show how the final formula 
follows from Theorem 0] applied to X given by (j!4j) with H replaced by 
H/y/M + N, d = 2, ai = lim M -»oc M/(M + iV) = 1/(1 + A), a 2 = A/(l + A); 
we only consider the case A > 1. From (JT5j) we have 

lim G hh */(m+n)(z) = ^-^tr a g(z) - = f(z), 

where Q{z) = diag(/(z), g(z)) satisfies equation (fT7|) with r/(diag(a, 6)) = 
diag(a 2 &, a\a); we used ()19|) to identify the answer with /. Thus 

2/ = 1 + Xzfg/(1 + A), zg = 1 + ^/(l + A). 

After a calculation which includes renormalization (dilation by 1/(1 + A)), 
this gives 

, /(^/(1 + A)) ^ + (l-A)-y/(z-l-A) 2 -4A 

^ - 1 + A - . 

which is the Cauchy transform of the Marchenko-Pastur law, see for example 
[22 page 102] □ 



5.2 Block Toeplitz matrices 

In this section we show how to use Theorem |2] to identify the limit laws of 
block Toeplitz matrices that were considered in ^2j- The interesting feature 
of this example is that Q is non-diagonal. In this example we present only de- 
tails for the case d = 3; in the figures we also show our result and comparison 
with simulation results for the cases d = 4 and d = 5. 

Suppose A, B, C are independent selfadjoint N x N matrices with i.i.d. 
entries of unit variance, and consider the block Toeplitz matrix: 

A B C 
BAB. 
C B A 
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Then as iV — > oo the spectral law of X converges almost surely to a deter- 
ministic law, see [12] • Since the limit does not depend on the distribution 
of the entries, we can assume that they are Gaussian and apply Theorem 
121 This determines the Cauchy transform of the limit law as the normalized 
trace of Q; the latter is determined by Equation (jl2j) . For Q of the form 



./ 





h 





9 





h 





./ 



1] acts as follows: 



„(5) = 3 



'2 / + g g + 2h 

2f+g+2h 
g+2h 2f+g 



Since the above pattern of Q is preserved under the mapping D I 3 + r](D) ■ 
D, we know that our wanted solution must have this form. So f!12l) gives the 
following system of equations: 

g (/ + h) + 2 {p + h 2 ) 



zf 



zg 



zh 



1 + 



g(g + 2(f + h)) 
3 

4fh + g (f + h) 



(22) 
(23) 
(24) 



This system of equations simplifies to two quadratic equations which can be 
solved by symbolic software; the density is a mixture of the semicircle law and 
another curve, but the exact solutions were not helpful in the selection of the 
appropriate root. Fig. |2(a)| illustrates the histogram and the theoretical den- 
sity obtained numerically from applying the Stieltjes inversion formula on the 
Cauchy transform G(z) = (2f(z) + g(z))/3, jj,{t) = — - lim s ^ (t + is). 
The numerical solution is calculated by solving the resulting system of equa- 
tions using Newton's method. The histograms are calculated for 100 matrices 
with Gaussian blocks of size 100 x 100. 

In Figures [2(b)1 and|2(c)| we compare the simulation and the result of our 
method for the cases d = 4 with 
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A 
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and d = 5 with 



A 


B 


C 


D 


E 


B 


A 


B 


C 


D 


C 


B 


A 


B 


C 


D 


C 


B 


A 


B 


E 


D 


C 


B 


A 



It might be of interest to point out that the double limit, first as block 
sizes iV — > oo and then d — > oo of block Toeplitz matrices gives the semicircle 
law. This can be seen by passing to d X d Toeplitz matrices with semicircular 
entries. Since the combinatorial estimates from (22] apply to this case, the 
moments are given as sums over pair partitions. For semicircular elements, 
only non-crossing partitions contribute, and each such partitions contributes 
1. 



5.3 Wishart matrix of rectangular block Toeplitz ma- 
trices 

In this section the introduced problem of calculating the spectrum for a 
MIMO channel is addressed, where it is of interest to calculate the distribu- 
tion of the Wishart matrix of rectangular block Toeplitz matrices. First for 
the case ut = n R , a proposition is presented that simplifies the process to 
calculate the desired spectrum. In the Example subsection, first the spec- 
trum for a MIMO system where ny = n R is calculated. Then the problem 
for the riR = 2nx is tackled and the results are presented. 

Assuming tit = Ur = N and denoting the channel matrix with Hn, the 
following proposition draws the way to calculate the eigenvalue distribution 
of H* N H N . 

5.3.1 Proposition 

As N — > oo, d = 2K + L — 1, n = Nd the spectral law of H^H* N /n converges 
weakly with probability one to the deterministic probability measure which 
is a mixture of K densities with Cauchy-Stieltjes transform 

1 - 

lim G HH * /n (z) = -V/^). 
N— >oo a 

i=i 

Functions fj are each a Cauchy transform of a probability measure and the 
following conditions hold. 



29 



1. f 3 = f K+1 .j for 1 < j < K 

2. There exist Cauchy transforms gi,g 2 , ■ ■ ■ ,gx+L-i of some probability 
measures such that gj = gx+L-j for 1 < j < K + L — 1 with the 
property that 

G(z) = diag(/i, f 2 , ■ ■ ■ , f K , 9i, 92,---, gK+L-i) 
satisfies (fT7|) with 

' lHk=i 9j+k, ttl<j<K, 

\ Etf h if K + 1 < j < K + mm{K, L}, 

3 ££f-K-L + i /* if ^ + ^ < J < and L < tf, 

3 Ef=i /* if 2if < j < K + L and L > 

I \ Etj-K-L+i fh if ^ + max{^, L} < j < 2K + L - 1. 



»7(0) 



3J 



Proof. We omit the proof of almost sure convergence, see Section 15.41 To 
verify that Q(z) is diagonal and its diagonal entries satisfy the symmetry 
condition we compute the diagonal entries of rj(D). For 1 < j < K we have 

L-l 



fc=0 

For 1 < j < mm{K, L} we have 

1 j 

[^( L ')]^+i,A-+i = j^2 d k,k- 
k=l 

For min{K, L} < j < max{i^, L} we have 



[rj{-L>)\K+j,K+j - \ 1 K J 

dLi=i d M L > K. 



For max{i^, L} < j < K + L — 1 we have 

K 



1 * 



d 

k=j-L+l 

We note that the symmetry conditions fj = fx+i-j and gj = gx+L-j are 
preserved under the mapping D I d -\- r](D) ■ D. Since r\ maps diagonal 
matrices into diagonal matrices, our solution Q must be diagonal and the 
diagonal entries must satisfy the symmetry conditions as claimed. □ 
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5.3.2 Examples 



In this part, first the developed theorems are used toward calculation of 
spectrum of a MIMO system with ISI (L = 4) and frame size of 4 (K = 4) 
while riT = ur. Then, for the simpler case of K = 2 and L = 2 when 
tir = 2rix (or = 2n#) the spectrum is calculated. 
Consider 



H 



N 



A B 
A 













c 

B 
A 




D 

C 
B 
A 




D 

C 
B 





D 

C 






D 



where A, B, C, D are independent non-self adjoint Gaussian iV x iV-random 
matrices. It is also assumed that the impulse response of the channel from 
any transmitter antenna to any receiver antenna is identical and equal to 
[ 1 1 1 1 ] . In this case K = L = 4, d = 11, 

Q{z) = diag(/i, / 2 , / 2 , fl, gi, 0}, 03, 04, 03, 02, 0l), 

and (j!7|) yields the following system of equations 



Zfl = 1 + 2/l(01 +02 + 03 +04)/H 

zf 2 = l + zf 2 (g 2 + 2g 3 + g A )/ll 

zgi = 1 + 2/15-1/11 

zg 2 = l + z(/i + / 2 WU 

^ 3 = l + 2(/i + 2/ 2 )53/H 

254 = l + 2z(/ 1 + / 2 )0 4 /H. 



(25) 
(26) 
(27) 
(28) 
(29) 
(30) 



The limiting Cauchy transform is Ghh*{z) = (/1 + /2V2. After eliminating 
from the equations we get 



1 1 

1 1 

J2 + IwT 



1 



1 



1 



11-/1 
+ 



/a 
2 



11 -A -2/2 11-2A-2/2 



+ 



/ 2 11-/!- 2/ 2 11 - 2/i - 2/ 2 



Mathematica failed to solve this system in a reasonable time. Fortunately, 
the Newton's algorithm can be applied to the equations (|23|l - (|3*n|) and its 
answer matches the simulations (Fig. 0J). 
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On the other hand, the channel matrix for a MIMO system with n# = 2ht 
and frame length of K = 2 over a channel with ISI (L = 2) is as follows: 



H 



N 



ABO 
A B 



where A and B are independent non-self adjoint 2N x N Gaussian random 
matrices. While one may calculate the spectrum of HH* using Theorem 
it is also possible to consider 
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N 



A 


B 





C 


D 








A 


B 





C 


D 



(31) 



for A, B, C and D independent non-self adjoint N x N Gaussian random 
matrices and calculate the spectrum of HH* using Theorem El The spectrum 
of HH* is calculated using the Newton's method which is depicted in Fig. 0J 



5.4 General remarks on convergence 

In Theorems 121 and 01 the convergence means weak convergence of the ex- 
pected values of the empirical laws of the eigenvalues, which is sufficient to 
identify the limit. In applications, the appropriate mode of convergence is 
weak convergence of random empirical measures with probability one. Such 
convergence means that for a "typical" realization of a random matrix, the 
limiting law as determined from Theorems El and H] is a good approximation 
when the size of the random matrix is large enough. 

Results of this type are readily available for symmetric independent blocks 
with i.i.d. entries, see [12]. In the Gaussian case one can easily extend 
such results to non-Hermitian and dependent square blocks by noting the 
following: 

• If A consists of i.i.d. Gaussian random variables, then A = A' / \/2 + 
A"/(V2i) where A' = (A + A*)/y/2 and A" = i{A - A*)/y/2 are in- 
dependent, Hermitian, and their off-diagonal entries are i.i.d Gaussian 
with the same variance as the entries of A. 

• If A, B are jointly Gaussian with respective entries that share the same 
correlation a (and are otherwise i.i.d), then B = A + yl — a 2 B', where 
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B 1 = (A — aB)/y/l — a 2 is independent of A and has i.i.d. entries of 
the same variance as the entries of A. 

For a Hermitian block matrix X consisting of square blocks, these two trans- 
formations allow us to express tr n (X fc ) as a linear combination of the polyno- 
mials in independent Gaussian and Hermitian matrices. Under appropriate 
normalization (i.e. when the matrices are constructed from an infinite i.i.d 
family of unit variance, and the matrices are scaled by the square root of the 
dimension) , the latter converge almost surely, see for example j23] . Since the 
limiting law is compactly supported, weak convergence of empirical measures 
holds with probability one. 

5.5 General remarks on the numerics 

The main inherent problem is that we have to find the "right" root. As we do 
not have a general automated recipe for doing this, we choose our solution by 
inspection. Theorem |21 guarantees uniqueness of an analytic solution under 
condition ([ll|h but the validity of this condition is difficult to judge for a 
numerical solutions. It seems that in general there is only one solution which 
corresponds to a valid probability measure. Furthermore, as it is quite easily 
seen, the fact that Q is essentially a generating power series in operator- 
valued moments implies that each diagonal entry of Q is itself the Cauchy 
transform of some probability measure. Thus the solution we are looking for 
must exhibit the very strong property that all diagonal entries of G{z) yield 
a valid probability distribution when z approaches the real line. 

Our experience with this is reflected quite adequately in the following 
remark of Edelman and Rao [25J on the same kind of problem with their 
random matrix calculator: "Irrespective of whether the encoded probability 
measure is compactly supported or not, the — 1/z behaviour of the real part 
of Stieltjes transform (the principal value) as z — > ±00 helps identify the 
correct root. In our experience, while multiple root curves might exhibit this 
behavior, invariably only one root will have an imaginary branch that, when 
normalized, will correspond to a valid probability measure. Why this always 
appears to be the case for the operational laws described is a bit of a mystery 
to us." 
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6 Conclusion 



We use operator-valued free probability theory and develop a method to 
calculate the limiting spectra for block matrices. Using this method we 
compute the eigenvalue distribution for block-matrices arising in wireless 
communications. The agreement with simulation results is excellent. We 
expect our method to have a much wider applicability. 
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channel impulse response 2] 2j 2] 2J 
between TX1 andRX2: [ h 1 h 2 h } ■■ ■ h L ] 



Figure 1: Block diagram of a MIMO system with ISI. 
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(c) 5 x 5 



Figure 2: Results for block Toeplitz matrices from Section 5.2. The his- 
togram in each figure stands for the simulation result (100 realization of 
matrices with 100 x 100 Gaussian blocks) while the darker curve represents 
the numerical results of the proposed method. 
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Figure 3: Superimposed theoretical density of the eigenvalues of complex 
normal H n H^/n for a channel with ISI L = 4 and a MIMO system tir = tit 
with frame length of K = 4 over its histogram for iV = 100, based on 100 
realizations. 




0.5 1 1.5 



Figure 4: Superimposed theoretical density of the eigenvalues of complex 
normal H n H*/n for a channel with ISI L = 2 and a MIMO system tir = 2ut 
with frame length of K = 2 over its histogram for iV = 100, based on 100 
realizations. 
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